home *** CD-ROM | disk | FTP | other *** search
/ Eagles Nest BBS 8 / Eagles_Nest_Mac_Collection_Disc_8.TOAST / Developer Tools⁄Additions / MacScheme20 / Mathlib / pseries.scm < prev    next >
Encoding:
Text File  |  1989-04-27  |  6.1 KB  |  196 lines  |  [TEXT/????]
  1. ;;; $Header: pseries.scm,v 1.4 88/01/26 07:44:32 GMT gjs Exp $
  2. ;;;; PS.SCM: Power-series arithmetic using infinite streams.
  3.  
  4. (if-mit
  5.  (declare (usual-integrations = + - * /
  6.                  zero? 1+ -1+
  7.                  ;; truncate round floor ceiling
  8.                  sqrt exp log sin cos)))
  9.  
  10. ;;; Power series may be represented by (infinite) streams of
  11. ;;;  coefficients.
  12.  
  13. ;;; Thus we may represent the series
  14. ;;;         a0 + a1*x + a2*x^2 + a3*x^3 + ...
  15. ;;;  by the stream {a0 a1 a2 a3 ...}
  16.  
  17. ;;; In our series the coefficients are manipulated with the
  18. ;;;  procedures ADD, MUL, SUB, DIV, rather than +, *, -, /, so
  19. ;;;  that we may use rational arithmetic, complex arithmetic,
  20. ;;;  or polynomial arithmetic... as desired.
  21. ;;;  Here we specify the map:
  22.  
  23. (define add +)
  24. (define sub -)
  25. (define mul *)
  26. (define div /)
  27. (define raise expt)
  28. (define zero 0)
  29.  
  30. ;;; The following procedures provide a set of capabilities for
  31. ;;;  manipulating series.
  32.  
  33. ;;; c*(a0 + a1*x + a2*x^2 + a3*x^3 + ...)
  34. ;;;  = c*a0 + c*a1*x + c*a2*x^2 + c*a3*x^3 + ...)
  35.  
  36. (define scale-series
  37.   (lambda (c s)
  38.     (map-stream (lambda (x) (mul c x))
  39.         s)))
  40.  
  41. ;;; (a0 + a1*x + a2*x^2 + ...) + (b0 + b1*x + b2*x^2 + ...)
  42. ;;;   = (a0+b0) + (a1+b1)*x + (a2+b2)*x^2 + ...
  43.  
  44. (define add-series
  45.   (lambda (s1 s2)
  46.     (map-streams add s1 s2)))
  47.  
  48. (define sub-series
  49.   (lambda (s1 s2)
  50.     (map-streams sub s1 s2)))
  51.  
  52.  
  53. ;;; (a0 + a1*x + a2*x^2 + ...) * (b0 + b1*x + b2*x^2 + ...)
  54. ;;;   = a0*b0 + (a0*b1+a1*b0)*x + (a0*b2+a1*b1+a2*b0)*x^2 + ...
  55. ;;; Each coefficient of the result is formed by reversing an initial
  56. ;;;  segment of one series, multiplying it by the coefficients of an
  57. ;;;  initial segment of the other series, and accumulating the
  58. ;;;  products.
  59.  
  60. (define mul-series
  61.   (lambda (s1 s2)
  62.     (map-streams inner-product-segments
  63.              (map-stream reverse (initial-segments-series s1))
  64.          (initial-segments-series s2))))
  65.  
  66. (define initial-segments-series
  67.   (let ()
  68.     (define i-s-helper
  69.       (lambda (seg series)
  70.     (let ((nseg (cons (head series) seg)))
  71.       (cons-stream nseg
  72.                (i-s-helper nseg
  73.                    (tail series))))))
  74.     (lambda (s)
  75.       (i-s-helper '() s))))
  76.  
  77. (define inner-product-segments
  78.   (lambda (seg1 seg2)
  79.     (reduce add (map mul seg1 seg2))))
  80.  
  81. (define div-series
  82.   (let ()
  83.     (define (div-helper v0)
  84.       (lambda (un ws rvs)
  85.     (div (sub un (inner-product-segments ws rvs)) v0)))
  86.     (lambda (u v)
  87.       (define v0 (head v))
  88.       (define w
  89.     (cons-stream
  90.           (div (head u) v0)
  91.           (map-streams (div-helper v0)
  92.                        (tail u)
  93.                        (initial-segments-series w)
  94.                        (map-stream reverse
  95.                                    (initial-segments-series (tail v))))))
  96.       w)))
  97.  
  98.  
  99. (define (expt-series v a)        ;A a number.
  100.   (cond ((zero? (head v))
  101.          (error "Cannot raise series to power -- no A0 term"))
  102.         ((= (head v) 1)
  103.          (expt-monic-series v a))
  104.         (else 
  105.          (let ((v0^a (raise (head v) a)))
  106.            (map-stream (lambda (wi) (mul wi v0^a))
  107.                        (expt-monic-series (map-stream (lambda (vi) 
  108.                                                         (div vi (head v)))
  109.                                                       v)
  110.                                           a))))))
  111.  
  112. (define (expt-monic-series v a)
  113.   (define w
  114.     (cons-stream 1
  115.                  (map-streams (lambda (w s1 s2)  ;weighted inner product
  116.                                 (reduce add (map mul w (map mul s1 s2))))
  117.                               (expt-weights (+ a 1))
  118.                               (initial-segments-series (tail v))
  119.                               (map-stream reverse
  120.                                           (initial-segments-series w)))))
  121.   w)
  122.  
  123. (define (expt-weights a+1)
  124.   (define (sloop n)
  125.     (let ((a+1/n (/ a+1 n)))
  126.       (define (wloop k l)
  127.         (if (> k n)
  128.             l
  129.             (wloop (+ k 1) (cons (- (* k a+1/n) 1) l))))
  130.       (cons-stream (wloop 1 '()) (sloop (+ n 1)))))
  131.   (sloop 1))
  132.  
  133. (define deriv-series
  134.   (let ()
  135.     (define (deriv-iter s n)
  136.       (if (null? s)
  137.       '()
  138.       (cons-stream (mul n (head s))
  139.                (deriv-iter (tail s) (+ n 1)))))
  140.     (lambda (s) (deriv-iter (tail s) 1))))
  141.  
  142.  
  143. ;;; The integral of a series
  144. ;;;           a0 + a1*x + a2*x^2 + a3*x^3 + ...
  145. ;;;  is       c + a0*x + a1*x^2/2 + a2*x^3/3 + ...
  146. ;;;  and is returned by the procedure INTEGRATE-SERIES which
  147. ;;;  takes the "initial condition" c as a required argument. 
  148. ;;;  For technical reasons, we are unable to use INTEGRATE-SERIES 
  149. ;;;  with mutual-recursion as in
  150. ;;;    (define cos-series (integrate-series sin-series 1)) ;DOESN'T
  151. ;;;    (define sin-series (integrate-series cos-series 0)) ;WORK!
  152. ;;;  However, we can achieve the desired effect by postponing the
  153. ;;;  attachment of the constant term, as follows. We use the procedure
  154. ;;;  INTEGRAL-SERIES-TAIL which returns the indefinite integral
  155. ;;;  part of the integrated series, i.e., {a0 a1/2 a2/3 a3/4 ...}.
  156. ;;;  Now, the mutual-recursion above can be made to work:
  157. ;;;    (define cos-series (cons-stream 1 (integral-series-tail sin-series)))
  158. ;;;    (define sin-series (cons-stream 0 (integral-series-tail cos-series)))
  159.  
  160. (define (integrate-series series constant-term)
  161.   (cons-stream constant-term (integral-series-tail series)))
  162.  
  163. (define integral-series-tail
  164.   (let ()
  165.     (define integrate-helper
  166.       (lambda (s n)
  167.     (cons-stream (div (head s) n)
  168.              (integrate-helper (tail s) (+ n 1)))))
  169.     (lambda (series)
  170.       (integrate-helper series 1))))
  171.  
  172.  
  173.  
  174. ;;; A power series may be evaluated at a particular
  175. ;;;  point.  Given a stream that represents a series, the
  176. ;;;  following procedure will produce a procedure that 
  177. ;;;  when given an x will produce a stream of pairs -- the 
  178. ;;;  first element of each pair is a partial sum, and the next 
  179. ;;;  element of the pair is the next (first ignored) term in
  180. ;;;  the sum.  Note that this is a sequence, not a series.
  181.  
  182. (define partial-sums
  183.   (let ((partial-sums-helper
  184.      (lambda (x)
  185.        (define loop
  186.          (lambda (series sum xn)
  187.                (let ((term (mul xn (head series))))
  188.              (cons-stream (list sum term)
  189.                   (loop (tail series)
  190.                     (add sum term)
  191.                     (mul xn x))))))
  192.        loop)))
  193.     (lambda (series)
  194.       (lambda (x)
  195.         ((partial-sums-helper x) series 0 1)))))
  196.